suppressPackageStartupMessages({
library(readr)
library(tidyverse)
library(reshape)
library(popbio)
})
## Conflicts with tidy packages ----------------------------------------------
Katsuwonus pelamis (Linnaeus, 1758), also known as Skipjack Tuna, is a fish from the Order Perciformes, under the Scombridae family. This family includes all mackerels, tunas, and bonitos. As other species in this family, K. pelamis has a pelagic affinity, often found in aggregations between 0 and 260 m deep and is mainly distributed in tropical waters –where they spawn– but is common to subtropical and temperate waters. Scombrids are also characterized for their fast growth rates and high mobility.
Unlike other Scombrids, K. pelamis is relatively small, with a mean length of 80 cm (max up to 110 cm). The maximum reported weight has been 34.5 Kg, and the maximum reported age has been 12 years. It reaches the sexual maturity around one year and 43 cm.
It feeds mainly on fish and crustaceans, but squids and some mollusks also comprise part of their diet. Trophic ecology research indicates that they have a trophic level of TL = 4.4 ± 0.5 (M ± SE; Froese & Pauly, 2016). This species is commercially important and sustain a high fishing effort around the world.
Besides its populations being classified as Least Concern by CITES, there is a high uncertainty is population estimations (which may mask the real status) and some regions are already showing the first signs of overfishing for this species.
The species has been identified as Least Concern by CITES. FAO reports increased landing in the last years, likely due to an increase in effort.
Catches of skipjack tuna have been steadily increasing since 1950, reaching a global peak in 1991 at 1 674 970 t. (FAO, 2017, http://www.fao.org/fishery/species/2494/en ). On the Atlantic Ocean the catches peaked on 2013 at 255,729.78 t.
Nominal catches data for Skipjack Tuna on the Atlantic Ocean, reported by countries to the International Commision to the Conservation of the Atlantic Tuna (ICCAT). The data was updated on November 2016 and incorporates catches from 1950 to 2015.
The data includes a total of 3,711 records of catch data, specified for 66 years (1950 - 2015). The data also includes information by Fleet (149), country (66 Flags), and Party (39). ICCAT Area code is included to allow spatial identification. In order to control for different gear types, a column for Gear group is included. Stock source (Eastern and Western Atlantic) is also included.
After comments on the first assignment, we looked for data on effort for Skipjack fisheries on the Atlantic Ocean. We found data from ICCAT on the number of fishing hours associated with catch from years 2006 to 2015. We then explored the relationship between effort (fishing hours) and time (Fig. 1), by fitting a linear model where coefficients were estimated by Ordinary Least Squares with heteroskedastic-robust standard errors, and tested the significance of the slope and intercept coefficients. Neither the intercept or slope showed significant change through time (Table 1), indicating that fishing effort had remained relatively constant from 2006 to 2015 (Fig 1).
Fig and Reg Table
We then used a larger dataset that included catches from 1950 to 2015, by country (Fig 2), and calculated the total catches by year (Fig. 3). We counted the unique countries that participated in the fishery or for which there is data (Figure 4), and used this as a proxy of total fishing effort. We then normalized total catches by dividing them by the number of countries participating in the fishery each year (Fig. 5). Here, we can see an increase in catches even after normalizing for our proxy of effort. For the Eastern stock (ATE), we see a steady increase since 1950, resulting in an order of magnitude higher values for 2015, as compared to 1950. On the other hand, the Western stock (WTE) does not show an increase of similar magnitude. Interannual variability in the data is in the order of 5,000 tonnes for ATE, and 1500 for ATW. Tables 2 and 3 present the regression coefficients for CPUE as a function of time, for each stock.
Figs and Reg tables
Our data is composed by length data, which we transformed into expected age using the transformed Von Bertalanffy function \[ A = \frac {1}{-K} \times log(1-(\frac {L}{L_{inf}}))+t_0 \], where: A = age, K = growth rate, L = length, \(L_{inf}\) = asymptotic length. Therefore, our model use a age structure.
Mortality: \(Z = 1.69\) (Garbin & Castello, 2014)
Fecundity: \(a = -1.33354\), \(b = 3.238\) (from FISHBASE)
von Bertanalnffy growth parameters: \(L_inf = 102.0\), \(K = 0.55\), \(t_0 = -0.02\) (Uchiyama & Struhsaker, 1981)
DiagrammeR::grViz("
digraph boxes_and_circles{
# Define nodes
node [shape = circle
penwidth = 2,
fontsize = 24]
egg; a1; a2; a3; a4; a5; a6; a7; a8; a9; a10; a11; a12; a13
# Add edge statements
#Advaance classes
egg -> a1[label = s0]
a1 -> a2[label = s1]
a2 -> a3[label = s2]
a3 -> a4[label = s3]
a4 -> a5[label = s4]
a5 -> a6[label = s5]
a6 -> a7[label = s6]
a7 -> a8[label = s7]
a8 -> a9[label = s8]
a9 -> a10[label = s9]
a10 -> a11[label = s10]
a11 -> a12[label = s11]
a12 -> a13[label = s12]
# Fecundities
a1 -> egg[label = f1]
a2 -> egg[label = f2]
a3 -> egg[label = f3]
a4 -> egg[label = f4]
a5 -> egg[label = f5]
a6 -> egg[label = f6]
a7 -> egg[label = f7]
a8 -> egg[label = f8]
a9 -> egg[label = f9]
a10 -> egg[label = f10]
a11 -> egg[label = f11]
a12 -> egg[label = f12]
a13 -> egg[label = f13]
}
", height = 1000)
Figure 1 - Life-cycle diagram for Skipjack tuna (Katsuwonus pelamis). S indicates survival, and f indicates fecundity.
The model we have defined axplicitly models eggs as a cohort (age = 0). We will conssider deleting this in the future, and explicitly calculate recruits.
#Build a matrix with text
A <- matrix(0, 14, 14) #initial empty matrix with all 0
# By using the $$, we can embed latex equations that will then be rendered to loog better
# Populate matrix with mortality
for (i in 2:14){
A[i,i-1] <- paste0("$s_{", i-1, ",", i-2, "}$")
}
ages <- seq(0:13)-1
A[1,] <- paste0("$f_{", seq(1:14)-1,"}$")
knitr::kable(A, col.names = paste0("$a_{",ages,"}$"), row.names = F, caption = "Table 1- Symbolical population matrix A. Subindices represent their rows and columns, respectively. The inferior diagonal represents survivals, while the first row represents facundities.")
| \(a_{0}\) | \(a_{1}\) | \(a_{2}\) | \(a_{3}\) | \(a_{4}\) | \(a_{5}\) | \(a_{6}\) | \(a_{7}\) | \(a_{8}\) | \(a_{9}\) | \(a_{10}\) | \(a_{11}\) | \(a_{12}\) | \(a_{13}\) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| \(f_{0}\) | \(f_{1}\) | \(f_{2}\) | \(f_{3}\) | \(f_{4}\) | \(f_{5}\) | \(f_{6}\) | \(f_{7}\) | \(f_{8}\) | \(f_{9}\) | \(f_{10}\) | \(f_{11}\) | \(f_{12}\) | \(f_{13}\) |
| \(s_{1,0}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | \(s_{2,1}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | \(s_{3,2}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | \(s_{4,3}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 | \(s_{5,4}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 | 0 | \(s_{6,5}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 | 0 | 0 | \(s_{7,6}\) | 0 | 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | \(s_{8,7}\) | 0 | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | \(s_{9,8}\) | 0 | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | \(s_{10,9}\) | 0 | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | \(s_{11,10}\) | 0 | 0 | 0 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | \(s_{12,11}\) | 0 | 0 |
| 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | \(s_{13,12}\) | 0 |
# Define Von Bert parameters
# Garbin & Castello, 2014
# L_inf <- mean(c(80, 80, 87.12, 94, 97.9, 97.3, 112.34, 89.38, 92.5))
# K <- mean(c(0.32, 0.6, 0.22, 0.38, 0.14, 0.25, 0.14, 0.38, 0.16))
# Or from Us and Stiurskjgfas, 1981
L_inf <- 102
K <- 0.55
t_0 <- -0.02
# Just to be clear, we use L-inf and to from Ushisomething, and K from the review
# Define fecundity parameters
# From fishbase http://www.fishbase.se/Reproduction/FecundityList.php?ID=107&GenusName=Katsuwonus&SpeciesName=pelamis&fc=416&StockCode=121
fec_a <- -1.33354
fec_b <- 3.238
# Define mortality
m <- 0.63
z <- 1.39
# Convert length to age using von bertalanffy model, solving for t
length2age <- function(length, l_inf, K, t_o){
age <- (1/-K)*(log(1-(length/L_inf))) + t_o
return(age)
}
# Convert age to length using von bertalanffy model
age2length <- function(age, l_inf, K, t_o){
length <- l_inf*(1-exp(-K*(age-t_o)))
return(length)
}
#Convert length to fecundity (number of eggs)
fecundity <- function(length, a, b){
f <- 10^(a+(b*log10(length*10)))
return(f)
}
A <- matrix(0, 14, 14) #initial empty matrix with all 0
# Populate matrix with mortality
for (i in 2:14){
A[i,i-1] <- exp(-1-z)
}
# Populate matrix with fecundity
ages <- seq(0:13)-1
lengths <- age2length(ages, L_inf, K, t_0)
A[1,] <- fecundity(lengths, fec_a, fec_b)
A[is.na(A)] <- 0
A[2,1] <- 0.0000001
A[1,1] <- 0
colnames(A) <- ages
rownames(A) <- ages
knitr::kable(A, col.names = paste0("$a_{",ages,"}$"), row.names = F, caption = "Table 2- Population matrix A. The inferior diagonal represents survivals, while the first row represents facundities.")
| \(a_{0}\) | \(a_{1}\) | \(a_{2}\) | \(a_{3}\) | \(a_{4}\) | \(a_{5}\) | \(a_{6}\) | \(a_{7}\) | \(a_{8}\) | \(a_{9}\) | \(a_{10}\) | \(a_{11}\) | \(a_{12}\) | \(a_{13}\) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0e+00 | 1.657256e+07 | 7.026737e+07 | 1.294406e+08 | 1.758241e+08 | 2.072322e+08 | 2.27012e+08 | 2.38998e+08 | 2.461086e+08 | 2.502768e+08 | 2.527038e+08 | 2.541114e+08 | 2.549259e+08 | 255396708 |
| 1e-07 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 9.162970e-02 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 9.162970e-02 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 9.162970e-02 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 9.162970e-02 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 9.162970e-02 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 9.16297e-02 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 9.16297e-02 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 9.162970e-02 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 9.162970e-02 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 9.162970e-02 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 9.162970e-02 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.00000e+00 | 0.00000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 9.162970e-02 | 0 |
load("size_dist_1980.RData") #This is the n-vector for 1980, loads it to a vector called n
n <- n %>% {
.$N}%>%
c(rep(0, times = 8))
n_0 <- sum(A[1, 2:14]*n, na.rm = T)
n <- c(n_0, n)
n
## [1] 2.562588e+19 8.674537e+06 3.136387e+11 2.711014e+10 4.391505e+08
## [6] 3.880114e+06 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [11] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
project <- popbio::pop.projection(A, n, 15)
pop <- project$stage.vectors %>%
as.data.frame() %>%
mutate(Age = as.factor(seq(0:13)-1)) %>%
gather(Year, N, -Age) %>%
mutate(Year = as.numeric(as.character(Year))) %>%
select(Year, Age, N)
## Warning in mutate_impl(.data, dots): '.Random.seed' is not an integer
## vector but of type 'NULL', so ignored
ggplot(pop, aes(x = Year, y = log(N), color = Age)) +
geom_line() +
theme_bw()
Figure 2 - Population size through time, represented by ages.
\(\lambda =\) 1.4729341
plot(project$pop.changes, type = "b", xlab = "Iterations", ylab = "Lambda")
Figure 3 - Convergence of \(\lambda\) through time.
pop %>%
filter(Year == 14) %>%
ggplot(aes(x = Age, y = log(N))) +
geom_bar(stat = "identity") +
theme_bw()
Figure 4 - Stabe age structure, reached after 15 iterations in the model.
sensitivity(A) %>%
image2(box.offset=.1)
Figure 5 - Sensitivity matrix.
elasticity(A) %>%
image2(box.offset = 0.1)
Figure 6 - Elasticity matrix.
The \(\lambda = 1.47\) shows that the population is projected to grow at rate of 47%. The result of the sensitivity analysis suggest that is most effective to intervene on the early ages. Since the species mature early (between ages 1 and 2), there is no necessity of conserving the old/big individuals that have higher fecundity, because the younger ones (in quantity) can provide enough eggs to maintain the population.
The model assumes the mortality is constant through ages. Better estimates of mortality/survival for different ages would show the real population structure.
The model assumes that the population is density-independent. However, we believe that fecundity, recruitment and survival of this pelagic species is density dependent. Perhaps a first approach wil be to use a Beverton-Holt approach to calculate recruitment, although this assumes density dependence in all the model.
Currently, all organisms are harvested in the model (excluding eggs). An interesting approach would be to set a size limit that reduces total effort on small sizes to natural effort by excluding fishing effort. This could allow us to evaluate the effect of a minimum catch size, management intervention aligned with the observed elasticities.
Environmental variation is another aspect to consider. It is known that recruitment of fish can be affected by extreme environmental conditions, like El Nino. Further analysis could incorporate the effects of climate variability on recruitment.
Mortality: \(Z = 1.39\) and \(M = 0.63\) (Garbin & Castello, 2014)
Fecundity: \(a = -1.33354\), \(b = 3.238\) (from FISHBASE)
von Bertanalnffy growth parameters: \(L_{inf} = 102.0\), \(K = 0.55\), \(t_0 = -0.02\) (Uchiyama & Struhsaker, 1981)
# Define Von Bert parameters
# Garbin & Castello, 2014
# L_inf <- mean(c(80, 80, 87.12, 94, 97.9, 97.3, 112.34, 89.38, 92.5))
# K <- mean(c(0.32, 0.6, 0.22, 0.38, 0.14, 0.25, 0.14, 0.38, 0.16))
# Or from Us and Stiurskjgfas, 1981
L_inf <- 102
K <- 0.55
t_0 <- -0.02
# Just to be clear, we use L-inf and to from Ushisomething, and K from the review
# Define fecundity parameters
# From fishbase http://www.fishbase.se/Reproduction/FecundityList.php?ID=107&GenusName=Katsuwonus&SpeciesName=pelamis&fc=416&StockCode=121
fec_a <- -1.33354
fec_b <- 3.238
# Define mortality
m <- 0.63
z <- 1.39
# Convert length to age using von bertalanffy model, solving for t
length2age <- function(length, l_inf, K, t_o){
age <- (1/-K)*(log(1-(length/L_inf))) + t_o
return(age)
}
# Convert age to length using von bertalanffy model
age2length <- function(age, l_inf, K, t_o){
length <- l_inf*(1-exp(-K*(age-t_o)))
return(length)
}
#Convert length to fecundity (number of eggs)
fecundity <- function(length, a, b){
f <- 10^(a+(b*log10(length*10)))
return(f)
}
A <- matrix(0, 14, 14) #initial empty matrix with all 0
# Populate matrix with mortality
for (i in 2:14){
A[i,i-1] <- exp(-z)
}
# Populate matrix with fecundity
ages <- seq(0:13)-1
lengths <- age2length(ages, L_inf, K, t_0)
A[1,] <- fecundity(lengths, fec_a, fec_b)
A[is.na(A)] <- 0
A[2,1] <- 0.0000001
A[1,1] <- 0
colnames(A) <- ages
rownames(A) <- ages
knitr::kable(A, col.names = paste0("$a_{",ages,"}$"), row.names = F,
caption = "Table I - Population matrix A. The inferior diagonal represents survivals, while the first row represents facundities.")
| \(a_{0}\) | \(a_{1}\) | \(a_{2}\) | \(a_{3}\) | \(a_{4}\) | \(a_{5}\) | \(a_{6}\) | \(a_{7}\) | \(a_{8}\) | \(a_{9}\) | \(a_{10}\) | \(a_{11}\) | \(a_{12}\) | \(a_{13}\) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0e+00 | 1.657256e+07 | 7.026737e+07 | 1.294406e+08 | 1.758241e+08 | 2.072322e+08 | 2.270120e+08 | 2.389980e+08 | 2.461086e+08 | 2.502768e+08 | 2.527038e+08 | 2.541114e+08 | 2.549259e+08 | 255396708 |
| 1e-07 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0 |
load("size_dist_1980.RData") #This is the n-vector for 1980, loads it to a vector called n
n <- n %>% {
.$N}%>%
c(rep(0, times = 8))
n_0 <- sum(A[1, 2:14]*n, na.rm = T)
n <- c(n_0, n)
n
## [1] 2.562588e+19 8.674537e+06 3.136387e+11 2.711014e+10 4.391505e+08
## [6] 3.880114e+06 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
## [11] 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00
We projected the population, stratified by ages, for 30 years in the future (Figure 1). The population seems to have a positive trend, which shows the population is increasing. Note how one-year olds achieve a population size of log(N) = 40 by year 22 (dashed lines). We will use this reference point when comparing with the management intervention.
project <- popbio::pop.projection(A, n, 30)
pop <- project$stage.vectors %>%
as.data.frame() %>%
mutate(Age = as.factor(seq(0:13)-1)) %>%
gather(Year, N, -Age) %>%
mutate(Year = as.numeric(as.character(Year))) %>%
select(Year, Age, N)
ggplot(pop, aes(x = Year, y = log(N), color = Age)) +
geom_line() +
theme_bw() +
geom_hline(yintercept = 40, color = "black", linetype = "dashed", size = 1) +
geom_vline(xintercept = 22, color = "black", linetype = "dashed", size = 1)
Figure 1 - Population size through time, represented by ages.
This time we projected the population for 30 years in the future with a management intervention on the minimum size of capture (Figure 2). Our projection assumes total compliance on the intervention, which means the fishing pressure is equal to 0. Therefore, the total mortality (Z) is equal to the natural mortality (M) for individuals of age 1. By protecting individuals of age 1 (the age at which they reach maturity), we observe a faster recovery in the fishery. This can be observed by comparing the log(N) = 40 at year 22, which is now reached at year 18 for one-year olds. Additionally, in the context of the entire stock, Figure 3 shows that with the management intervention total skipjack population recovers faster.
A_min <- A
A_min[3,2] <- exp(-m)
knitr::kable(A_min, col.names = paste0("$a_{",ages,"}$"), row.names = F, caption = "Table I - Population matrix A, modifying mortality of age 1 organisms. The inferior diagonal represents survivals, while the first row represents facundities.")
| \(a_{0}\) | \(a_{1}\) | \(a_{2}\) | \(a_{3}\) | \(a_{4}\) | \(a_{5}\) | \(a_{6}\) | \(a_{7}\) | \(a_{8}\) | \(a_{9}\) | \(a_{10}\) | \(a_{11}\) | \(a_{12}\) | \(a_{13}\) |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0e+00 | 1.657256e+07 | 7.026737e+07 | 1.294406e+08 | 1.758241e+08 | 2.072322e+08 | 2.270120e+08 | 2.389980e+08 | 2.461086e+08 | 2.502768e+08 | 2.527038e+08 | 2.541114e+08 | 2.549259e+08 | 255396708 |
| 1e-07 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 5.325918e-01 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0.000000e+00 | 0 |
| 0e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 0.000000e+00 | 2.490753e-01 | 0 |
project_min <- popbio::pop.projection(A_min, n, 30)
pop <- project_min$stage.vectors %>%
as.data.frame() %>%
mutate(Age = as.factor(seq(0:13)-1)) %>%
gather(Year, N, -Age) %>%
mutate(Year = as.numeric(as.character(Year))) %>%
select(Year, Age, N)
ggplot(pop, aes(x = Year, y = log(N), color = Age)) +
geom_line() +
theme_bw() +
geom_hline(yintercept = 40, color = "black", linetype = "dashed", size = 1) +
geom_vline(xintercept = 22, color = "black", linetype = "dashed", size = 1) +
geom_vline(xintercept = 18, color = "red", linetype = "dashed", size = 1)
Figure 2 - Population size through time, represented by ages.
tot1 <- project$pop.sizes
tot2 <- project_min$pop.sizes
time <- seq(1:30)
data.frame(time, BAU = tot1, Int = tot2) %>%
gather(scenario, popsize, -time) %>%
ggplot(aes(x = time, y = log(popsize), color = scenario)) +
geom_point(size = 2) +
geom_line(size = 1) +
theme_bw() +
scale_color_brewer(palette = "Set1")
Figure 3 - Total population size through time.